from EOS_collection import CUBIC_EOS
import numpy as np

R = 8.314

def T_BOIL_CAL(Tc, Pc, w, Z, PS):

    #温度初始化
    T0 = 10
    T = 1000
    T1 = 0

    #求PR方程的参数
    PR = CUBIC_EOS(Pc, Tc, Z, w)
    b = PR.PR_b()

    #初始化
    VSV = 0
    VSL = 0

    #设定初始猜测温度的下限
    while abs(VSL - VSV) < 1e-5:

        #初始化迭代值
        a = PR.PR_a(T0)

        #计算气相的体积根
        VSV = R * T0 / PS
        VSVnew = PR.PR_V(PS, VSV, T0, a, b)
        while abs(VSVnew - VSV) / VSV > 1e-10:
            VSV = VSVnew
            VSVnew = PR.PR_V(PS, VSV, T0, a, b)

        #计算液相的体积根
        VSL = 2 * b
        VSLnew = PR.PR_V(PS, VSL, T0, a, b)
        while abs(VSLnew - VSL) / VSL > 1e-10:
            VSL = VSLnew
            VSLnew = PR.PR_V(PS, VSL, T0, a, b)

        #计算压缩因子
        ZSV = PS * VSV / R / T0
        ZSL = PS * VSL / R / T0

        #计算气相的逸度系数
        phiSV = PR.PR_PHI(PS, VSV, T0, a, b, ZSV)

        #计算液相的逸度系数
        phiSL = PR.PR_PHI(PS, VSL, T0, a, b, ZSL)

        #计算函数值
        F0 = phiSV - phiSL

        #缩小温度范围
        T0 = T0 + 1

    #重新初始化
    VSV = 0
    VSL = 0

    #设定初始猜测温度的上限
    while abs(VSL - VSV) < 1e-5:

        #初始化迭代值
        a = PR.PR_a(T)

        #计算气相的体积根
        VSV = R * T / PS
        VSVnew = PR.PR_V(PS, VSV, T, a, b)
        while abs(VSVnew - VSV) / VSV > 1e-10:
            VSV = VSVnew
            VSVnew = PR.PR_V(PS, VSV, T, a, b)

        #计算液相的体积根
        VSL = 2 * b
        VSLnew = PR.PR_V(PS, VSL, T, a, b)
        while abs(VSLnew - VSL) / VSL > 1e-10:
            VSL = VSLnew
            VSLnew = PR.PR_V(PS, VSL, T, a, b)

        #计算压缩因子
        ZSV = PS * VSV / R / T
        ZSL = PS * VSL / R / T

        #计算气相的逸度系数
        phiSV = PR.PR_PHI(PS, VSV, T0, a, b, ZSV)

        #计算液相的逸度系数
        phiSL = PR.PR_PHI(PS, VSL, T0, a, b, ZSL)

        #缩小温度范围
        T = T - 1

    #初始化逸度
    phiSV = -1
    phiSL = -2

    while abs(phiSV - phiSL) > 1e-5:

        a = PR.PR_a(T)

        #计算气相的体积根
        VSV = R * T / PS
        VSVnew = PR.PR_V(PS, VSV, T, a, b)
        while abs(VSVnew - VSV) / VSV > 1e-10:
            VSV = VSVnew
            VSVnew = PR.PR_V(PS, VSV, T, a, b)

        #计算液相的体积根
        VSL = 2 * b
        VSLnew = PR.PR_V(PS, VSL, T, a, b)
        while abs(VSLnew - VSL) / VSL > 1e-10:
            VSL = VSLnew
            VSLnew = PR.PR_V(PS, VSL, T, a, b)

        #计算压缩因子
        ZSV = PS * VSV / R / T
        ZSL = PS * VSL / R / T

        #计算气相的逸度系数
        phiSV = PR.PR_PHI(PS, VSV, T, a, b, ZSV)

        #计算液相的逸度系数
        phiSL = PR.PR_PHI(PS, VSL, T, a, b, ZSL)

        F = phiSV - phiSL

        #割线法迭代
        T1 = T - F * ( T - T0 ) / ( F - F0 )
        T0 = T
        T = T1
        F0 = F

    return T
